# sample code for ordinal logistic regression

nall <- 900
m<-4
run<-0
times<-1000
thetamean<-matrix(0,(m+3),1)
sdmean<-matrix(0,(m+3),1)
HL_o<-matrix(0,times,3)
AUC_o<-rep(0,times)
sc<-rep(0,times)
while(run<times)
{	
	
	
# generating data with sample size 1800	
	v1<-rnorm(nall)
	v2<-rnorm(nall)
	v3<-rbinom(nall,1,0.5)
	v4<-rbinom(nall,1,0.5)
	

	 
	 
	 
	xall<-matrix(1,nall,m)
	xall<-cbind(v1,v2,v3,v4)
	
	
	intcept<-matrix(1,nall,1)
	L1<--1*intcept+1*xall[,1]+1*xall[,2]+1*xall[,3]+1*xall[,4]
	L2<-0*intcept+1*xall[,1]+1*xall[,2]+1*xall[,3]+1*xall[,4]
	L3<-1*intcept+1*xall[,1]+1*xall[,2]+1*xall[,3]+1*xall[,4]
	p01<-exp(L1)/(1+exp(L1))
	p02<-exp(L2)/(1+exp(L2))
	p03<-exp(L3)/(1+exp(L3))
	
	



	yall<-matrix(0,nall,4)

	for(i in 1:nall)
		yall[i,] <- rmultinom(1,1,c(p01[i],(p02[i]-p01[i]),(p03[i]-p02[i]),(1-p03[i])))

	x<-xall[1:450,]
	y<-yall[1:450,]
	xtest<-xall[451:900,]
	ytest<-yall[451:900,]


# estimating on ordinal model using 900 sample

    n<-450
	theta0<-matrix(-1,(m+3),1)
	theta1<-matrix(0,(m+3),1)
	theta1[1,1]<-0
	theta1[2,1]<-0.1
	theta1[3,1]<-0.2
	
	
	k_cen<-0
	epsilon<-10^(-6) 
	
	while(max(abs(theta1-theta0))>epsilon&k_cen<20)
	{ 
	  theta0<-theta1
	  alpha1<-theta0[1,1]
	  alpha2<-theta0[2,1]
	  alpha3<-theta0[3,1]
	  beta<-theta0[seq(4,(m+3)),1]
	  p1<-exp(alpha1+x%*%beta)/(1+exp(alpha1+x%*%beta))
	  p2<-exp(alpha2+x%*%beta)/(1+exp(alpha2+x%*%beta))
	  p3<-exp(alpha3+x%*%beta)/(1+exp(alpha3+x%*%beta))
	  p1sq<-exp(alpha1+x%*%beta)/((1+exp(alpha1+x%*%beta))^2)
	  p2sq<-exp(alpha2+x%*%beta)/((1+exp(alpha2+x%*%beta))^2)
	  p3sq<-exp(alpha3+x%*%beta)/((1+exp(alpha3+x%*%beta))^2)
	  a1a2<-exp(alpha1)/(exp(alpha2)-exp(alpha1))
	  a2a1<-exp(alpha2)/(exp(alpha2)-exp(alpha1))
	  a2a3<-exp(alpha2)/(exp(alpha3)-exp(alpha2))
	  a3a2<-exp(alpha3)/(exp(alpha3)-exp(alpha2))
	  
	  # derivative for alpha1
	  da1<-sum(y[,1]*(1-p1)+y[,2]*(-a1a2-p1))
	  # derivative for alpha2
	  da2<-sum(y[,2]*(a2a1-p2)+y[,3]*(-a2a3-p2))
	  # derivative for alpha3
	  da3<-sum(y[,3]*(a3a2-p3)+y[,4]*(-p3))
	  # derivative for beta
	  db<-t(x)%*%(y[,1]*(1-p1)+y[,2]*(1-p2-p1)+y[,3]*(1-p3-p2)+y[,4]*(-p3))
	  # second derivative for alpha1
	  d2a1<-sum(y[,1]*(-p1sq)+y[,2]*(-a1a2*a2a1-p1sq))
	  # second derivative for alpha2
	  d2a2<-sum(y[,2]*(-a1a2*a2a1-p2sq)+y[,3]*(-a2a3*a3a2-p2sq))
	  # second derivative for alpha3
	  d2a3<-sum(y[,3]*(-a2a3*a3a2-p3sq)+y[,4]*(-p3sq))
	  # second derivative for alpha1 and alpha2
	  d2a12<-sum(y[,2]*(a1a2*a2a1))
	  # second derivative for alpha2 and alpha3
	  d2a23<-sum(y[,3]*(a2a3*a3a2))
	  # second derivative for alpha1 and beta
	  d2a1b<-t(x)%*%((y[,1]+y[,2])*(-p1sq))
	  # second derivative for alpha2 and beta
	  d2a2b<-t(x)%*%((y[,2]+y[,3])*(-p2sq))
	  # second derivative for alpha3 and beta
	  d2a3b<-t(x)%*%((y[,3]+y[,4])*(-p3sq))
	  # second derivative for beta
	  temp<-y[,1]*(-p1sq)+y[,2]*(-p2sq-p1sq)+y[,3]*(-p3sq-p2sq)+y[,4]*(-p3sq)
	  d2b<-t(x)%*%diag(c(temp))%*%x
	  # first derivative for theta 
	  temp1<-matrix(0,3,1)
	  temp1[1,1]<-da1
	  temp1[2,1]<-da2
	  temp1[3,1]<-da3
	  dall<-rbind(temp1,db)
	  # second derivative for theta
	  temp2<-matrix(0,3,3)
	  temp2[1,1]<-d2a1
	  temp2[2,2]<-d2a2
	  temp2[3,3]<-d2a3
	  temp2[1,2]<-d2a12
	  temp2[2,1]<-d2a12
	  temp2[2,3]<-d2a23
	  temp2[3,2]<-d2a23
	  temp3<-cbind(d2a1b,d2a2b,d2a3b)
	  temp2<-rbind(temp2,temp3)
	  temp3<-rbind(t(temp3),d2b)
	  d2all<-cbind(temp2,temp3)
	  d2all<-(-d2all)
	  
 	  # newton-raphson updata
	  theta1<-theta0+solve(d2all+diag(0.000001,(m+3)))%*%(dall)
	  k_cen<-k_cen+1
	  cat("\n\nk_cen=",k_cen)
	}
	if (k_cen<20)
	{	
	    run<-run+1
		cat("\n\nrun=",run)
		
		thetamean<-thetamean+theta1
			
		hattheta<-theta1
		alpha1<-hattheta[1,1]
		alpha2<-hattheta[2,1]
		alpha3<-hattheta[3,1]
		beta<-hattheta[seq(4,(m+3)),1]
		p1<-exp(alpha1+x%*%beta)/(1+exp(alpha1+x%*%beta))
		p2<-exp(alpha2+x%*%beta)/(1+exp(alpha2+x%*%beta))
		p3<-exp(alpha3+x%*%beta)/(1+exp(alpha3+x%*%beta))
		p1sq<-exp(alpha1+x%*%beta)/((1+exp(alpha1+x%*%beta))^2)
		p2sq<-exp(alpha2+x%*%beta)/((1+exp(alpha2+x%*%beta))^2)
		p3sq<-exp(alpha3+x%*%beta)/((1+exp(alpha3+x%*%beta))^2)
		a1a2<-exp(alpha1)/(exp(alpha2)-exp(alpha1))
		a2a1<-exp(alpha2)/(exp(alpha2)-exp(alpha1))
		a2a3<-exp(alpha2)/(exp(alpha3)-exp(alpha2))
		a3a2<-exp(alpha3)/(exp(alpha3)-exp(alpha2))
		
		
		
		#######CORRECT WAY FOR SD#################
		# second derivative for alpha1
		d2a1<-sum(y[,1]*(-p1sq)+y[,2]*(-a1a2*a2a1-p1sq))
		# second derivative for alpha2
		d2a2<-sum(y[,2]*(-a1a2*a2a1-p2sq)+y[,3]*(-a2a3*a3a2-p2sq))
		# second derivative for alpha3
		d2a3<-sum(y[,3]*(-a2a3*a3a2-p3sq)+y[,4]*(-p3sq))
		# second derivative for alpha1 and alpha2
		d2a12<-sum(y[,2]*(a1a2*a2a1))
		# second derivative for alpha2 and alpha3
		d2a23<-sum(y[,3]*(a2a3*a3a2))
		# second derivative for alpha1 and beta
		d2a1b<-t(x)%*%((y[,1]+y[,2])*(-p1sq))
		# second derivative for alpha2 and beta
		d2a2b<-t(x)%*%((y[,2]+y[,3])*(-p2sq))
		# second derivative for alpha3 and beta
		d2a3b<-t(x)%*%((y[,3]+y[,4])*(-p3sq))
		# second derivative for beta
		temp<-y[,1]*(-p1sq)+y[,2]*(-p2sq-p1sq)+y[,3]*(-p3sq-p2sq)+y[,4]*(-p3sq)
		d2b<-t(x)%*%diag(c(temp))%*%x
		#################################################
		
		
		
		
		#  second derivative for theta
		temp2<-matrix(0,3,3)
		temp2[1,1]<-d2a1
		temp2[2,2]<-d2a2
		temp2[3,3]<-d2a3
		temp2[1,2]<-d2a12
		temp2[2,1]<-d2a12
		temp2[2,3]<-d2a23
		temp2[3,2]<-d2a23
		temp3<-cbind(d2a1b,d2a2b,d2a3b)
		temp2<-rbind(temp2,temp3)
		temp3<-rbind(t(temp3),d2b)
		d2all<-cbind(temp2,temp3)
		d2all<-(-d2all)
		# variance covariance matrix
		cov_matri<-solve(d2all)
		# standard error
		sd<-sqrt(diag(cov_matri))
		sdmean<-sdmean+sd
		
		
		# Hosmer-Lemeshow test for ordinal estimation 
		ntest<-450
		p1<-exp(alpha1+xtest%*%beta)/(1+exp(alpha1+xtest%*%beta))
	    p2<-exp(alpha2+xtest%*%beta)/(1+exp(alpha2+xtest%*%beta))
	    p3<-exp(alpha3+xtest%*%beta)/(1+exp(alpha3+xtest%*%beta))
		
		ytest1<-ytest[,1]
		ytest1<-matrix(ytest1,ntest,1)
		ytest2<-ytest[,1]+ytest[,2]
		ytest2<-matrix(ytest2,ntest,1)
		ytest3<-ytest[,1]+ytest[,2]+ytest[,3]
		ytest3<-matrix(ytest3,ntest,1)
		
		sub<-as.integer((ntest)/10)
		n_tot<-matrix(sub,10,1)
		n_tot[10,1]<-(ntest)-9*sub
		
		#test1
		hltable<-cbind(ytest1,p1)
		hltable<-hltable[order(hltable[,2]),]
		n_obs<-matrix(0,10,1)
		n_est<-matrix(0,10,1)
		
		for(i in 1:9)
		{
			temp<-hltable[(1+(i-1)*sub):(i*sub),1]
			n_obs[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),2]
			n_est[i]<-sum(temp)
			
		}
		temp<-hltable[(1+9*sub):(ntest),1]
		n_obs[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),2]
		n_est[10]<-sum(temp)
		
		HL1<-sum((n_obs[n_est!=0&n_est!=n_tot]-n_est[n_est!=0&n_est!=n_tot])^2/(n_est[n_est!=0&n_est!=n_tot]*(1-n_est[n_est!=0&n_est!=n_tot]/n_tot[n_est!=0&n_est!=n_tot])))
		HL_o[run,1]<-HL1
		#test2
		hltable<-cbind(ytest2,p2)
		hltable<-hltable[order(hltable[,2]),]
		n_obs<-matrix(0,10,1)
		n_est<-matrix(0,10,1)
		
		for(i in 1:9)
		{
			temp<-hltable[(1+(i-1)*sub):(i*sub),1]
			n_obs[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),2]
			n_est[i]<-sum(temp)
			
		}
		temp<-hltable[(1+9*sub):(ntest),1]
		n_obs[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),2]
		n_est[10]<-sum(temp)
		
		HL2<-sum((n_obs[n_est!=0&n_est!=n_tot]-n_est[n_est!=0&n_est!=n_tot])^2/(n_est[n_est!=0&n_est!=n_tot]*(1-n_est[n_est!=0&n_est!=n_tot]/n_tot[n_est!=0&n_est!=n_tot])))
		HL_o[run,2]<-HL2
		#test3
		hltable<-cbind(ytest3,p3)
		hltable<-hltable[order(hltable[,2]),]
		n_obs<-matrix(0,10,1)
		n_est<-matrix(0,10,1)
		
		for(i in 1:9)
		{
			temp<-hltable[(1+(i-1)*sub):(i*sub),1]
			n_obs[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),2]
			n_est[i]<-sum(temp)
			
		}
		temp<-hltable[(1+9*sub):(ntest),1]
		n_obs[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),2]
		n_est[10]<-sum(temp)
		
		HL3<-sum((n_obs[n_est!=0&n_est!=n_tot]-n_est[n_est!=0&n_est!=n_tot])^2/(n_est[n_est!=0&n_est!=n_tot]*(1-n_est[n_est!=0&n_est!=n_tot]/n_tot[n_est!=0&n_est!=n_tot])))
		HL_o[run,3]<-HL3
		
		
		
		
		
		
		# AUC score for ordinal estimation
		p1<-exp(alpha1+xtest%*%beta)/(1+exp(alpha1+xtest%*%beta))
	    p2<-exp(alpha2+xtest%*%beta)/(1+exp(alpha2+xtest%*%beta))
	    p3<-exp(alpha3+xtest%*%beta)/(1+exp(alpha3+xtest%*%beta))
		
		ytest1<-ytest[,1]
		ytest1<-matrix(ytest1,ntest,1)
		ytest2<-ytest[,1]+ytest[,2]
		ytest2<-matrix(ytest2,ntest,1)
		ytest3<-ytest[,1]+ytest[,2]+ytest[,3]
		ytest3<-matrix(ytest3,ntest,1)
		
		
		
		# AUC1
		auctable<-cbind(ytest1,p1)
		tot<-auctable
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc1<-sumrank/(length(pone)*length(pzero))
		
		# AUC2
		auctable<-cbind(ytest2,p2)
		tot<-auctable
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc2<-sumrank/(length(pone)*length(pzero))
		
		# AUC3
		auctable<-cbind(ytest3,p3)
		tot<-auctable
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc3<-sumrank/(length(pone)*length(pzero))
		
		AUC<-(auc1+auc2+auc3)/3
		AUC_o[run]<-AUC
		
		
        		
		# score test for proportional odds assumption
		temp<-theta1[c(1,seq(4,(m+3)),2,seq(4,(m+3)),3,seq(4,(m+3))),1]
		theta0<-matrix(temp,3*(m+1),1)
		beta1<-theta0[seq(1,(m+1)),1]
		beta2<-theta0[seq((m+2),(2*m+2)),1]
		beta3<-theta0[seq((2*m+3),(3*m+3)),1]
		# add intercept column to old x
		intcept<-matrix(1,n,1)
		x<-cbind(intcept,x)
		p1<-exp(x%*%beta1)/(1+exp(x%*%beta1))
		p2<-exp(x%*%beta2)/(1+exp(x%*%beta2))
		p3<-exp(x%*%beta3)/(1+exp(x%*%beta3))
		p1sq<-exp(x%*%beta1)/((1+exp(x%*%beta1))^2)
		p2sq<-exp(x%*%beta2)/((1+exp(x%*%beta2))^2)
		p3sq<-exp(x%*%beta3)/((1+exp(x%*%beta3))^2)
		b1b2<-exp(x%*%beta1)/(exp(x%*%beta2)-exp(x%*%beta1))
		b2b1<-exp(x%*%beta2)/(exp(x%*%beta2)-exp(x%*%beta1))
		b2b3<-exp(x%*%beta2)/(exp(x%*%beta3)-exp(x%*%beta2))
		b3b2<-exp(x%*%beta3)/(exp(x%*%beta3)-exp(x%*%beta2))

		# derivative for beta1
		db1<-t(x)%*%(y[,1]*(1-p1)+y[,2]*(-b1b2-p1))
		# derivative for beta1
		db2<-t(x)%*%(y[,2]*(b2b1-p2)+y[,3]*(-b2b3-p2))
		# derivative for beta1
		db3<-t(x)%*%(y[,3]*(b3b2-p3)+y[,4]*(-p3))
		# Expectation of second derivative for beta1
		temp<-(p1)*(-p1sq)+(p2-p1)*(-b1b2*b2b1-p1sq)
		d2b1<-t(x)%*%diag(c(temp))%*%x
		# Expectation of second derivative for beta2
		temp<-(p2-p1)*(-b1b2*b2b1-p2sq)+(p3-p2)*(-b2b3*b3b2-p2sq)
		d2b2<-t(x)%*%diag(c(temp))%*%x
		# Expectation of second derivative for beta3
		temp<-(p3-p2)*(-b2b3*b3b2-p3sq)+(1-p1)*(-p3sq)
		d2b3<-t(x)%*%diag(c(temp))%*%x
		# Expectation of second derivative for beta1 and beta2
		temp<-(p2-p1)*(b1b2*b2b1)
		d2b12<-t(x)%*%diag(c(temp))%*%x
		# Expectation of second derivative for beta2 and beta3
		temp<-(p3-p2)*(b2b3*b3b2)
		d2b23<-t(x)%*%diag(c(temp))%*%x
		# first derivative for theta
		dall<-rbind(db1,db2,db3)
		# Expectation of second derivative for theta
		zero<-matrix(0,(m+1),(m+1))
		temp2<-rbind(d2b1,d2b12,zero)
		temp3<-rbind(d2b12,d2b2,d2b23)
		temp4<-rbind(zero,d2b23,d2b3)
		d2all<-cbind(temp2,temp3,temp4)
		d2all<-(-d2all)
		sc[run]<-t(dall)%*%solve(d2all)%*%dall
		
		
			
	}
}

